Dispersion Properties of Multilayered 
Metal-Dieletric Metamaterials 



Alexey A. Orlov^ Ivan V. lorsh^^, Pavel A. Belov\ Yuri S. Kivshar^^ 

^ National Research University of Information Technologies, Mechanics and Optics 
197101, Kronverksky pr. 49, St. Petersburg, Russian Federation 
^St. Petersburg Academic University - Nanotechnology Research and Education Centre, 
St. Petersburg 194021, Russia Federation 
^Nonlinear Physics Centre, Research School of Physics and Engineering, Australian National 
University, Canberra ACT 0200, Australia 

alexey. orlov @ phoi. ifmo. ru 

Abstract: We study dispersion properties of layered metal-dielectric me- 
dia having different layers thicknesses ratios. Plotting dispersion diagrams 
and isofrequency contours, we find that strong nonlocalality is inherent 
property of such periodic structures. We introduce different level of losses 
and analyze complex modes of the metamaterial demonstrating that it 
operates in a regime with infinite numbers of eigenmodes, with nonlocality 
slightly affected by losses. 
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1. Introduction 

Multilayered structures are well-known basic structures in optics. They have simple geometry, 
therefore such structures can be studied with ease. Nevertheless, this simplicity hasn't pre- 
vented multilayered structures from catching a new wave of attention last years. Optical meta- 
materials formed by multilayered metal-dielectric nanostructures (MMDN) were found to have 
striking electromagnetic properties. An impressive example is the fact that the structure, be- 
ing treated as an indefinite medium [l], is able to experience effective metric signature changes 
from Minkowski-like ( — h ++) to 2T physics ( h+) A number of applications exploit- 
ing unique properties of the multilayered metamaterial were proposed, namely subwavelength 
imaging |[3tiTT1l. nanolithography |[T2l . optical nanocircuits |[T3l , 3D negative refraction O, 
and invisibility flS^ . Recently, multilayered metamaterials have become widely considered as 



Fig. 1. Scheme of the structure. 



a reahzation of the so-called hyperbolic media |[T6l with resulting ultra-high values for the 
Purcell factor ifTTlfTSl . 

It is known that the structure under consideration demonstrates a strong nonlocal re- 
sponse 1 19 |. Previously, we have shown in 1201 that under certain conditions two extraordinary 
waves instead of one can be excited in the metamaterial by incident TM-polarized light beam. 
The existence of such an additional light wave is a clear manifestation of the presence of spa- 
tial dispersion in the metamaterial, and it validates strong nonlocal nature of the structure. In 
the recent paper [21] this phenomenon was treated in the framework of arising field of nonlo- 
cal transformation optics. An important question of will the nonlocality and the beam splitting 
phenomenon still remain in real structures with introduction of losses emerges though. Conse- 
quently, it is needed to study complex modes of the multilayered metamaterial and clarify an 
effect of losses. 

Local effective medium model along with transfer matrix method ll22ll are used in order to 
describe the structure. The former models the structure depicted in Fig. [T] with layers permittiv- 
ities £i , £2 and thicknesses di , J2 as the uniaxial anisotropic crystal of the following form: 

/£^ 0\ ^11 

Ceff = £,1 , 

\0 £j £^ 

In the case when £nand £± are of opposite sign, i.e. < 0, a material is called indefinite fl] 
or hyperbolic ||23l|24l. Since hyperbolic metamaterials allow high-^ modes they may have very 
high photonic density of states ||25]|26l . 

However, it should be noticed that the local effective medium model is not accurate and 
used in the present paper only to show nonlocal behaviour of the multilayered metamaterial by 
comparison between results provided by the local model and the actual description. The nonlo- 
cal effective medium model (see Ref. |27|) is capable of providing analytic results coinciding 
with ones of transfer matrix method which is the most correct electromagnetic description of 
multilayers. 

Analysis of MMDN was accomplished in the work |28 | for lossless metamaterials. However, 
proper attention was not devoted to the imprescriptible attribute of the structure - nonlocality. 
Inapplicability of quasi-static approximation, that is effective medium model, wasn't demon- 
strated as well. We will start below from a demonstration of the fact that the effective medium 
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Fig. 2. Dispersion diagrams for different layers configurations: (a) di < d2, (b) di = d2, 
and (c) di > d2. Convergence dynamics with the decrease of the period D is shown, (d) 
Fine structure of the dispersion curves of the metamaterial in the case of thin metal. Four 
modes are observed to co-exist in the same frequency range. 



model can't be applied for accurate description of the multilayered metal-dielectric metamate- 
rials even in the deep subwavelength limit. 

2. Dispersion Diagrams 

First of all, we study electromagnetic properties of the lossless structure by means of dispersion 
diagrams and isofrequency contours. Losses will be introduced in Sec. HI 

For MMDN with two layers forming its period, using effective medium model and transfer 
matrix method we obtain two dispersion relations ()L)(k), where k = {kx^ky^Q) is the wave vector. 
TM polarization is assumed to take advantage of plasmonic behaviour. The first relation is given 
by effective medium model: 

e|l £± \ c J 

and another one is given by transfer matrix method, assuming that k^^ is the wave number in 
the /-th layer: 

co^{KD) = co5{ki'U,)cos{kPd2) - I + ^] sm{ki'U,)5m{kPd2). (3) 

In this section graphic representation of the dispersion dependencies is chosen in the form of 
dispersion diagrams for waves travelling along j-direction, i.e. along the layers with kx = and 
presented in Fig. [21 Isofrequency contours will be shown in Sec.O 



General view of the structure is presented in Fig.Oa. The structure consists of metal and di- 
electric layers and indefinite in all directions. Period of the structure is Z) = Ji + J2 = 62.5nm, 
where di and J2 are the layers thicknesses. Dielectric function for the dielectric layers is con- 
stant and equals to 4.6, while metal layers are described by dielectric function of the Drude 
form: 

2 

£2(0) = 1 - /"^^ (4) 

(o[(o-\-ir) 

with the plasma wavelength = 2nc/ (Op = AD and the damping coefficient P. The multilay- 
ered metamaterial with 3 different ratios J1/J2 of layers thicknesses (2/3, 1 and 3/2, respec- 
tively) and different values of D were chosen in order to demonstrate different behaviours of 
the structure. 

The case of the full period equal to 62.5 nm has been considered in Ref. |29|. Here we also 
consider structures with periods equal to Z)/2,Z)/4, and Z)/8 and compare their behaviour with 
one of the full period. Astonishingly, the behaviour observed for the full period that is 10 times 
smaller than the wavelength remains qualitatively the same for structures having 2, 4 and 8 
times smaller periods (i.e. 20, 40 and 80 times smaller than wavelength), as it is shown in Fig.O 
The convergence to an asymptote corresponding to surface plasmon-polariton (SPP) resonance 
becomes weaker as the branches reach asymptote at larger values of ky, but a topology of the 
branches remains the same. A pronounced tendency is such that the decrease of the period 
pushes one of dispersion curves closer to one predicted by the effective medium model for 
small ky. For large ky that convergence to effective medium model is weaker. Another dispersion 
branch tends to converge to a straight line corresponding to a volume plasmon with the decrease 
of the period. Volume plasmon is at the frequency where £|| = and it is marked out with 
green thin line in Fig.[2l Positions of these two lines, effective medium model line and volume 
plasmon line, depend on the /J2 ratio. That means that the convergence also depends also on 
the thicknesses ratio. In the case of equal thicknesses the convergence is faster than for different 
thicknesses. However, it is clear that even in the quasi-static limit correspondence between the 
effective medium model predictions and the actual dispersion properties of the metamaterial 
is not satisfactory, especially for high-/: waves with relatively large ky which are mainly used 
in sub wavelength imaging. We have also discovered from dispersion diagram analysis that the 
metamaterial in the case of thin metal {d\ > J2) possesses a fine structure shown in Fig.Od and 
supports four modes co-existing in the same frequency range. 

3. Isofrequency Contours 

Another powerful tool to analyze dispersion relation (o{k) is analysis of isofrequency contours, 
which are obtained if one fixes the frequency CO and plots kx{ky) dependence. Isofrequency con- 
tours for three different thicknesses ratios of multilayered metamaterials at certain frequencies 
are presented in Figs.lMSl You can find isofrequency contours for the whole optical frequency 
range in the video attachment of the present paper. 

For the configuration corresponding to the case of thin metal (di > J2), which has the most 
unusual dispersion properties, isofrequency contours are presented in Fig. [3] for 8 different fre- 
quencies. Examination of the contours reveals that at certain frequencies they have very unusual 
topology. For Z)/A = 0.066 (see Fig.[3la), two hyperbolic isofrequency contours predicted by 
effective model are well reproduced by MMDN. Deviations from the effective medium model 
appear only at large kx close to the edges of the Brillouin zone, where the homogenization pro- 
cedure is forbidden since the wavelength in the material becomes comparable to its period. For 
Z)/A = 0.087 (see Fig. Ob), where the frequency is closer to the SPP resonance, in addition to 
the hyperbolic contours MMDN demonstrates an elliptic-like contour right around the center 
of Brillouin zone. This region of the Brillouin zone is usually treated as one where the quasi- 
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Fig. 3. Isofrequency contours for the case of di > d2. Group velocities are shown for real- 
istic contours by arrows. 
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Fig. 4. Isofrequency contours for the case of di = d2. 
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Fig. 5. Isofrequency contours for the case of di < d2. 



static approximation is to be applicable since eigen waves have small wave numbers, i.e. small 
spatial variation. However, inside layers fields vary dramatically and this leads to emergence 
of nonlocality. At slightly higher frequency Z)/A = 0.089 (see Fig.[3lc), elliptical and hyper- 
bolic contours observed at lower frequencies are transformed into a pair of crossing ellipses. 
In this case, the effective medium model predicts completely flat isofrequency contour which 
is perfect for realization of canalization regime for sub wavelength imaging ||5l[8l, but it is not 
in agreement with realistic dispersion properties of MMDN. At the frequency D/X = 0.096 
(see Fig. [3]d), an elliptical isofrequency contour predicted by effective medium model is ac- 
companied by two hyperbolic contours, which disappear at higher frequencies, for example, at 
Z)/A = 0.106 (see Fig. [He). However, they appear once again at slightly higher frequencies, 
e.g. Z)/A = 0.106 (see Fig. [Hf). Presence of elliptic isofrequency contours is a clear sign of 
nonlocality in the MMDN. AiD/X = 0.108 (see Fig. Og) the topology of isofrequency con- 
tours is similar to ones at Z)/A = 0.089 (the crossing ellipses), but at slightly higher frequency 
D/X = 0.109 (see Fig. Oh) the contours acquire even more unusual, pitcher- like form, which 
also has nothing to do with effective medium model predictions. 

For the rest of thicknesses ratios when the thicknesses are equal or dielectric is thinner than 
metal isofrequency contours are presented in Figs. HI [5] There are not new topological structures 
that are absent in the case of thin metal, however, hyperbolic-like and elliptical contours are still 
there. One should note flat isofrequency contours in Fig. lit that appear at the frequency where 
= oo. It is known that flat isofrequency contours are in use for sub wavelength image canal- 
ization 1811301. Also, in the case of equal thicknessess plotted in Fig.|4lwe see that isofrequency 
contours of effective medium model are in a quite good agreement with realistic ones. 



4. Complex Band Structure 

Dispersion relation given by Eq. [3] allows us to study complex modes of the structure immedi- 
ately complex ky is considered. There are plenty of methods complex dispersion equation can 
be solved by, including a method based on the argument principle theorem 1311 . operator theory 
methods l32l . and numerical FEM methods I33II34I . However, in order to find roots of Eq.[3]in 
the complex plane we employ a combination of change of sign method and Newton's one. 

In Figs. O [7] dispersion diagrams of lossless (F 0) and lossy structures showing complex 
modes are presented. We show positive real part of the wave vector with corresponding imag- 
inary part. Due to the periodicity and unboundedness of the structure its complex eigenvalue 
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Fig. 6. Complex band structure of the structure with di < d2 in the absence of losses (a) and 
with the losses (b). Imaginary parts of the modes are on the left, real parts are on the right. 
Low-loss modes are marked with colors, gray curves corresponds to high degenerative 
modes. 



spectrum contains an infinite number of modes. However, only modes having low imaginary 
parts of its wavenumber are of interest since they can travel significant distances. For example, 
modes having Im{ky) < n/D are able to propagate to distances more than D/ln. We restrict 
imaginary part of ky by the value of 8 tt/Z) in our calculations. 

In the absence of losses the structure supports the set of waves having the wavevectors ±^ 
and =b^*. Thus, there are forward and backward waves along with their complex conjugations, 
and imaginary parts of the modes are symmetric in Figs. [6^, [7^, . It follows from the fact that if 
dielectric permittivities of layers are real then the right part of EqOis a holomorphic function 
whose restriction to the real numbers is real- valued. 

The most important three modes are marked in Figs. [6l [7] with different styles and colors 
while high-order complex modes are in gray color. Effective medium model results are pre- 
sented in the figures as well (thin curves). Modes showed in Fig.[2]of Sec.[2]can be recognized 
by their zero imaginary part in regions of propagation. Where in Fig. [2l there were band gaps, 
here we have fully complex modes. At the SPP resonance frequency real part of the propa- 
gating mode goes to infinity, while imaginary part degenerates at this point what can be seen 
on the example of mode II in the case of d\ < J2 and on the mode I in the case of d\ > J2. 
One can also observe an effect of decay switching in high-order modes: after a mode under- 



a) 0.13 

0.12 
0.11 
0.1 
0.09 
0.08 







' 1 ' E \. 






1 






1 i 






1 






III E 






1 
























\ / 












ni/ 11*^-^^— 




















Hi 






1 T 1 

1 ♦ - • 1 1 
• • 























-4 -2 2 4 

Im(ky)T>lTi 



2 4 

Re(kJD/Ti 



-6 -4 -2-0.3 0.3 2 4 6 

Im(kJD/Ti 



2 4 

Re(kJD/n 




Fig. 7. Complex band structure of the structure with di > d2 in the absence of losses (a) 
and with the losses (b). Imaginary parts of the modes are on the left, real parts are on the 
right. 



goes the SPP resonance its imaginary part jumps from one branch to another. In the case of 
di < d2 that is presented in Fig. [6^ there is a very special mode III that is fully imaginary in the 
presented frequency range, and its real part is dispersionless and equals to zero. This feature 
is very unique since that mode is not sensitive to the resonance. Its presence is highly related 
to effective medium model as the latter has the resonance at another frequency and at higher 
frequencies converges to the entirely imaginary mode III while at low frequencies it converges 
to the mode II. 

However, zero losses aren't realistic and unachievable in practice. Thus, we introduce at- 
tenuation in metal to study effect of losses in metal layers on electromagnetic properties of 
MMDN. Increasing losses smoothly we find that it changes step-like behaviour of the modes at 
the SPP resonance. That is, a specific form of these steps is determined by the value of losses in 
MMDN. With change of F modes go through a number of transformations ||35]| . Mechanism of 
modes transformation is exposed in Fig.[8l When losses tend to zero, at the resonance positive 
imaginary branches jump to their left neighbours with growth of the frequency, while negative 
branches jump to their right neighbours. Then, increase of the losses level forces branches to 
modify their topology and after a number of metamorphoses at the resonance both positive and 
negative imaginary branches jump to their right neighbours with growth of the frequency. 

This way we come to realistic level of losses. F = 1.734 • 10^^^"^ is supposed in our calcula- 
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tions. Since now we analyze complex modes in the case of realistic losses comparing complex 
band structure with the lossless case to show how losses affect the behaviour of the structure. 
In the presence of losses, as soon as the damping coefficient F in the dielectric function of 
metal £2 becomes positive non-zero and £2 becomes a complex number, the right part of Eq.[3] 
can be complex- valued with real ky. Consequently, the modes confluence is removed and the 
conjugated waves with are no longer supported by the structure. However, correspondence 
between positive and negative real and imaginary parts can be quite diverse because of the 
absence of symmetry in the imaginary parts of eigemodes (see Figs. [Tt, [lb). 

In Figs. [TJ), [6]3 we still have decay switching effect that is seen better in high degenerative 
modes. It's expected for modes we have observed in Figs. [7^, [6^: losses of a mode are grown 
up dramatically when it enters into a band gap, but for complex modes such behaviour is non- 
obvious. Also, as soon as one moves from the real axis to the complex plane every mode 
becomes complex one and it becomes hard to distinguish between propagating and evanescent 
modes. The case of di > J2 is shown in Fig.[6l3. Just that very case corresponds to simultaneous 
presence of forward and backward modes in the same frequency range (Z)/ A : 0. 105 — 0. 108) for 
propagation along the layers. That is, the beam splitting phenomenon revealed in ll2Ql remains 
even after introduction of realistic level of losses. Decay switching effects are observed there 
for high degenerative modes and the modes I, III while the mode II has narrow propagation 
domain with high losses in the rest of it. 

Finally, it's noteworthy that effective medium model predicts resonance behaviour for imagi- 
nary part of ky. In reality, however, the structure is functioning in a regime with infinite number 
of modes which imaginary parts have step-like behaviour at the resonance. 

5. Field Profiles 

In order to obtain deeper understanding of the nature of the studied complex eigenmodes of the 
system, we have plotted profiles of the magnetic and electric fields of these modes. Fig.|9]shows 
the profiles of the 3 branches for the case of the thicker metal at the frequency Z)/A = 0.0875. 
We first note that the magnetic field of the mode I averaged over each layer is equal to zero. 
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Fig. 9. Field profiles in the case of di < d2 for three modes at the frequency of D/X = 
0.0875. Solid curves correspond to the real part of the magnetic field, dashed curves corre- 
spond to the imaginary part of the electric field. 




Fig. 10. Field profiles of high-order modes corresponding to dispersion curves shown in 
Fig.md. 
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Fig. 11. Field profiles in the case of di > d2 for three modes at the frequency of D/X = 
0.0875. Solid curves correspond to the real part of the magnetic field, dashed curves corre- 
spond to the imaginary part of the electric field. 



It means that the mode I is effectively longitudinal. Such modes are well known in plasma 
physics as Langmuir modes 1361 . It has been also recently shown that such exotic modes may 
exist in the waveguides with metal cladding and core made of hyperbolic media 1371 . We can 
see however, that the longitudinal nature of the mode holds only within the effective media 
approximation and the local magnetic field inside each layer does not vanish. Modes II and III 
are conventional transverse modes. 

Another interesting feature that can be noted on the profiles is the complicated field distribu- 
tion of the electric field inside the metal layer for mode I. Conventionally the modes existing 
in metal-dielectric structures are coupled surface plasmon resonances at the individual metal- 
dielectric interfaces. Thus, there may exist no more than one extremum of the field distribution 
function in each layer. This condition however does not hold for the case of the complex modes. 
Indeed, if we consider the modes with real and imaginary part of the propagation constant, the 
expression for the transverse component of the wave vector in each layer reads: 

kx = yjeikl-k] = ^ Eikl - Kt{kyY ^lm{kyf - 2iKt{ky)lm{ky) . (5) 

It is clear that when lm(ky) is large enough, the k^ gains large real part and the mode starts to 
propagate inside the structure. This leads to nonzero values of the transverse component of the 
Poynting vector inside the layers and to the complex shape of the field distribution inside the 
layers. To illustrate how the imaginary part of ky affects the shape of the field inside the layers 
we have plotted the profiles of the high order complex modes [see Fig.[TO||. We can see that 
while we are switching to higher order modes, having the larger imaginary part of ky the field 
profile shape becomes more and more complicated. These modes can be regarded as a specific 
type of coupled waveguide modes. The waveguide mode condition can be roughly estimated as 
ki^xdi = Tin, where / = 1,2 and ^ - is an integer. Thus, higher switching between the complex 
modes can be regarder as switching between different waveguide modes in the structure. 



We have also plotted the profiles of the three mode branches for the thinner metal case [see 
Fig- El- We can see that in this case we also have one longitudinal mode I and two transverse 
modes. In this case however the mode II has the complicated structure of the electric field as 
the one defined by the larger values of Im(^_y). 

6. Conclusion 

Using dispersion diagrams along with isofrequency contours we have analyzed different char- 
acteristic configuration of MMDN both in the absence and in the presence of losses. Study of 
real- valued eigenmodes has revealed the impossibility of its description in terms of the local 
effective medium model even in the deep subwavelength limit. With the introduction of realis- 
tic losses, complex band structure of MMDN was revealed. Generally, nonlocality in MMDN 
is slightly affected by losses. Thus, in real structures with losses variety of spatial dispersion 
effects is expected. 
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